      PROGRAM NAD2(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT,PUNCH)
C                            PROGRAM NAD2
C
C NAD2 IS A FINITE DIFFERENCE PROGRAM EMPLOYING MURRAY LANDIS TRANSFORMATION
C FOR TWO PHASE BINARY SYSTEMS. THIS PROGRAM ALLOWS FOR PLANAR, CYLINDRICAL,
C OR SPHERICAL GEOMETRY. B STANDS FOR THE SOLUTE OR FILAMENT ELEMENT WHILE A
C STANDS FOR MATRIX ELEMENT. ALL CONCENTRATIONS ARE READ IN ATOM FRACTION B,
C ALL DISTANCES IN CENTIMETERS AND TIME IN SECONDS.
C
      INTEGER PCS
      DIMENSION XCAP(7), YCAP(7), TYTLE(8),FI(8), ABLE(16), XSN(501)
      DIMENSION RESULT(2),CSN(501  ),XN(501)
      COMMON/STORE/NX,CN(51 ),DN(51), ND,CP(501),DP(501) ,DMAX
      DATA RESULT/10HFID FOR TW,10HO PHASE   /
C
C PSEUDO AND FONTS ARE SYSTEM SUBROUTINES NECESSARY FOR PLOTTING.
C
      CALL PSEUDO
      CALL FONTS(1)
C
C                             FORMAT STATEMENTS
C
  100 FORMAT(2X,//)
  101 FORMAT (2I4)
  102 FORMAT(8A10)
  103 FORMAT(5E14.5)
  104 FORMAT ( 1H1 ,//)
  105 FORMAT ( 1H , 8A10)
  107 FORMAT(//2X, 36HGRID CHANGE WAS INITIATED WHEN TIME= ,E11.4,9H SEC
     1ONDS. /2X,55HFOLLOWING DATA PERTAINS TO SOLUTION BEFORE GRID CHANG
     2E. /2X,22HBETA-PHASE THICKNESS = ,E11.4,13H CENTIMETERS./2X,
     322HALPHA-PHASE THICKNESS= ,E11.4,13H CENTIMETERS. /)
  108 FORMAT(//2X,E12.4,5X,16HATOMIC MASS OF A/2X,E12.4,5X,16HATOMIC MAS
     1S OF B /2X,E12.4,5X,20HPAR. MOLAR VOL. OF A /2X,E12.4,5X, 20HPAR.
     2MOLAR VOL. OF B //)
  110 FORMAT(F8.5,2X,E14.8)
  111 FORMAT(2X,I4,2E15.6,2X,I4,2E15.6,2X,I4,2E15.6)
  112 FORMAT(4E15.6)
  113 FORMAT(4X,13HCONCENTRATION ,5X,21HDIFFUSION COEFFICIENT/(3X,F11.4,
     112X,E12.4))
  114 FORMAT(//5X,I4,13X,72HNXA-NUMBER OF X VALUES TO BE SET UP IN THE A
     1LPHA PHASE OF DIFFUSION ZONE /5X,I4,13X, 71HNXB-NUMBER OF X VALUES
     2 TO BE SET UP IN THE BETA PHASE OF DIFFUSION ZONE /2X,E15.6,5X,
     327HTHB-THICKNESS OF BETA PHASE/ 2X,E15.6,5X, 28HTHA-THICKNESS OF A
     4LPHA PHASE /2X,E15.6,5X, 25HTOTAL LENGTH OF DIF. ZONE/2X,E15.6,5X,
     530HFINAL LENGTH OF DIFFUSION ZONE /2X,E15.6,5X,50HINITIAL ESTIMATE
     6D DIFFUSION DISTANCE IN BETA PHASE /5X,F4.1,13X,3HPCS)
  115 FORMAT ( //3( 3X, 3HNO., 5X, 8HDISTANCE,  5X, 11HCOMPOSITION,1X))
  116 FORMAT(//1H ,5X,23HSTARTING DIFFUSION TIME , 7X,1H=,E15.8,1X,
     16HSECS =  ,E15.8, 5H MINS /6X,20HFINAL DIFFUSION TIME,
     210X, 1H=, E15.8,    7H SECS = ,E15.8, 5H MINS
     3/6X,18HINITIAL THICKNESS=, E11.4,17H FOR BETA-PHASE, ,
     4E11.4,18H FOR ALPHA-PHASE
     5/6X,18HFINAL THICKNESS  =, E11.4,17H FOR BETA-PHASE, ,
     6E11.4,18H FOR ALPHA-PHASE
     7/6X, 24HLENGTH OF DIFFUSION ZONE ,1H=,E15.8,2X/ 6X,
     842HTOTAL DIST. ALPHA/BETA INTERFACE HAS MOVED , 1X,1H=,E15.8 /
     96X, 25HINITIAL AREA UNDER CURVE= ,E15.8/6X, 25HFINAL AREA UNDER CU
     ARVE  = ,E15.8 /6X, 53HRATIO OF FINAL TO INITIAL AMOUNT OF B IN BET
     BA PHASE = ,E15.8//)
  124 FORMAT (4F10.5)
  132 FORMAT(//2X,27HINITIAL COMPOSITION PROFILE )
  133 FORMAT(//2X,25HPROFILE AFTER GRID CHANGE )
  209 FORMAT (F10.2)
  247 FORMAT(10F6.2)
  248 FORMAT(I2,2X,7A10)
  251 FORMAT(F6.2)
C
C            PLOT PARAMETERS
C
C XAL    LENGTH OF X-AXIS IN INCHES.
C YAL    LENGTH OF Y-AXIS IN INCHES.
C DXN    NUMBER OF INTERVALS ON X-AXIS WHERE THE TICK MARKS ARE TO BE PLACED
C DYN    NUMBER OF INTERVALS ON Y-AXIS WHERE THE TICK MARKS ARE TO BE PLACED
C XMIN   VALUE CORRESPONDING TO FIRST TICK MARK ON THE X-AXIS
C YMIN   VALUE CORRESPONDING TO FIRST TICK MARK ON THE Y-AXIS
C XIN    SPACING BETWEEN SUCCESSIVE TICK MARKS ON THE X-AXIS
C YIN    SPACING BETWEEN SUCCESSIVE TICK MARKS ON THE Y-AXIS
C XALM   MAX LENGTH OF X-AXIS IN INCHES.
C SPS    OPTION FOR SAVING PLOT SPECIFICATIONS.
C        0, AUTO SCALING BY INCREASING XIN SO THAT THE MAXIMUM
C           DISTANCE XMAX IN THE PROFILE IS LESS THAN XALM.
C        1, SAVE THE PLOT SPECIFICATIONS BY IGNORING THE POINTS
C           FOR WHICH X IS GREATER THAN THE LAST X-TICK MARK.
C NCX    NUMBER OF CHARACTERS IN THE X-CAPTION.
C XCAP   X-CAPTION
C NCY    NUMBER OF CHARACTERS IN THE Y-CAPTION.
C YCAP   Y-CAPTION
C TYTLE  TYTLE TO BE PRINTED IN THE PLOT.
C NA     PARAMETER DICTATING NUMBER OF SOLUTIONS PLOTTED ON EACH GRAPH.
C        1, ONE SOLUTION PER GRAPH
C        2, PLOT ALL SOLUTIONS ON THE SAME GRAPH.
C
      READ(5,247) XAL,YAL,DXN,DYN,XMIN,YMIN,XIN,YIN,XALM,SPS
      READ(5,248)NCX,XCAP
      READ(5,248) NCY,YCAP
      READ(5,102) TYTLE
      READ(5,101) NA
C
C ABLE      IDENTIFICATION FOR SOLUTION.
C AMW,BMW   THE ATOMIC MASSES FOR THE PURE SPECIES OF ALPHA AND BETA PHASES
C AMV,BMV   PARTIAL MOLAR VOLS. OF PURE SPECIES OF ALPHA AND BETA PHASES
C CZERO     CONCENTRATION IN THE MATRIX.
C CPRIME    CONCENTRATION IN THE B RICH REGION.
C           CPRIME,CZERO ARE USED TO SET UP A STARTING PROFILE
C NDA,NDB   NO. OF D VALUES IN THE ALPHA AND BETA PHASES
C CN        CONCENTRATION, ATOM FRACTION B.
C DN        DIFFUSION COEFFICIENT, IN CENTIMETER SQUARED PER SECOND.
C ND        NUMBER OF ( CONCENTRATION, DIFFUSION COEFFICIENT ) SETS IN TABLE.
C
C SPECIAL NOTE.. SOLUBILITY LIMITS FOR EACH PHASE ARE TAKEN
C AS THE END POINT CONCENTRATIONS (CN) FOR THAT PHASE.
C
C PCS       PARAMETER DICTATING FORM OF SOLUTION
C           0 , PLANAR GEOMETRY
C           1 , CYLINDRICAL GEOMETRY
C           2 , SPHERICAL GEOMETRY
C RC        INITIAL LENGTH OF TOTAL DIFFUSION ZONE.
C THB       INITIAL THICKNESS OF BETA.
C RCF       FINAL LENGTH OF TOTAL DIFFUSION ZONE.
C RD        ESTIMATED DIFFUSION ZONE FOR BETA.
C NXA,NXB   NO. OF X-LOCATIONS IN THE ALPHA AND BETA PHASES WHERE
C           THE CONC. IS TO BE EVALUATED.THE NO. OF X-INTERVALS IN
C           THE ALPHA AND BETA PHASES IS EQUAL TO NXA-1 AND NXB-1
C NSP       OPTION FOR READING STARTING CONCENTRATION PROFILE
C           1, READ AN INITIAL PROFILE
C           0, THE PROGRAM INITILISES THE PROFILE WITH CPRIME IN THE
C           FILAMENT, AND CZERO IN MATRIX.
C NO        PARAMETER DICTATING FORM OF OUTPUT.
C           1, WRITTEN ONLY
C           2, WRITTEN AND PUNCHED
C FI        VARIABLE FORMAT FOR PUNCHED OUTPUT OF CONCENTRATION PROFILE.
C
      READ (5,102)  ( ABLE(I), I=1,16 )
      READ (5,124)  AMW,BMW
      READ (5,124)    AMV,BMV
      READ (5,124) CZERO, CPRIME
      READ (5,101)  NDA,NDB
      ND=NDA+NDB
      READ(5,110) (CN(I),DN(I),I=1,ND)
      READ(5,101) PCS
      READ(5,112) RC,THB,RCF,RD
      READ(5,101)  NXA,NXB
      READ(5,101)  NSP
      READ(5,101)  NO
      READ(5,102)  (FI(I),I=1,8)
      WRITE(6,104)
      WRITE (6,105) ( ABLE(I), I=1,16 )
      WRITE(6,108) AMW,BMW,AMV,BMV
      WRITE(6,113) (CN(I),DN(I),I=1,ND)
C
C THIS SEGMENT INSPECTS THE TABLE OF DIFFUSION COEFFICIENTS TO SEE IF IT IS
C ARRANGED IN DESCENDING ORDER OF CONCENTRATION AND, IF NOT, ARRANGES IT SO.
C
      IF ( CN(1).GE.CN(ND) )  GO TO 4
      DO 2 I=1,ND
      CP(I)=CN(I)
    2 XSN(I) = DN(I)
      DO 3 I=1,ND
      N = ND+1-I
      CN(I)= CP(N)
    3 DN(I) = XSN(N)
    4 CONTINUE
C
C THE SOLUBILITY LIMITS ARE TAKEN AS FOLLOWS.
C
      NDBP1= NDB+1
      CBA= CN(NDB)
      DCBA= DN(NDB)
      CAB= CN(NDB+1)
      DCAB= DN(NDB+1)
      DMAXA=0.
      DMAXB=0.
      DO 1 I=1,NDB
    1 IF(DN(I) .GT.DMAXB)  DMAXB=DN(I)
      DO 7 I=NDBP1,ND
    7 IF(DN(I) .GT.DMAXA)  DMAXA=DN(I)
      DMAX= AMAX1(DMAXA, DMAXB)
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS IN ATOM FRACTION TO
C VOLUME FRACTION.
C
      DO 394 I=1,ND
      CN(I)=(CN(I)*BMV)/(AMV+CN(I)*(BMV-AMV))
  394 CONTINUE
      CZERO  =( CZERO  *BMV)/ (AMV+CZERO  *(BMV-AMV))
      CAB    =( CAB    *BMV)/ (AMV+CAB    *(BMV-AMV))
      CBA    =( CBA    *BMV)/ (AMV+CBA    *(BMV-AMV))
      CPRIME =( CPRIME *BMV)/ (AMV+CPRIME *(BMV-AMV))
      CDELTA=CBA-CAB
      T=0.
      IF (NSP.EQ.0)  GO TO 9
C
C THE FOLLOWING SEGMENT IS NOT USED UNLESS A STARTING PROFILE IS READ IN.
C ABLE      IDENTIFICATION OF PROFILE.
C FI        VARIABLE FORMAT FOR READING IN PROFILE.
C T         TIME OF PROFILE IN SECONDS.
C ICAW      OPTION TO READ THE PROFILE IN ATOM FRACTION OR WEIGHT FRACTION.
C           +1, WEIGHT FRACTION
C           -1, ATOM FRACTION
C NP        NUMBER OF ( X, CONCENTRATION ) SETS IN PROFILE.
C CSN       CONCENTRATION.
C XSN       DISTANCE, IN CENTIMETERS.
C
      READ (5,102)  ( ABLE(I), I=1,16 ),  ( FI(I), I=1,8 )
      READ (5,103)  T
      READ(5,101) ICAW
      READ (5,101)  NP
      READ (5,FI)   (CSN(I  ),XSN(I),I=1,NP)
      WRITE(6,FI) (CSN(I  ),XSN(I),I=1,NP)
      IF(ICAW.GT.0)GO TO 622
      DO 621 I=1,NP
  621   CP(I) =(CSN(I  )*BMV)/(AMV+CSN(I  )*(BMV-AMV))
      GO TO 623
  622 DO 620 I=1,NP
      CP(I)   =(CSN(I  )*(BMV/BMW))/((CSN(I  )*(BMV/BMW))+((1.0-CSN(I  )
     1)*(AMV/AMW)))
  620 CONTINUE
  623 CONTINUE
      CT=CDELTA/2.0
      KT=0
   29 KT=KT+1
      KT1=KT+1
      CTEST=CSN(KT)- CSN(KT1)
      IF (CTEST.LT.CT) GO TO 29
      THB=(XSN(KT)+XSN(KT1))/2.0
      WRITE (6,105)  ( ABLE(I), I=1,16 )
    9 CONTINUE
      TSTART=T
      TFINAL=T
      TSMIN=T/60.
      TFMIN= T/60.
      GEF= 4.* (1.+FLOAT(PCS))
C
C THE FOLLOWING SEGMENT SETS UP THE X ARRAY TO BE USED IN THE SOLUTION.
C NX IS THE NUMBER OF X-LOCATIONS WHERE THE COMPOSITION IS TO BE CALCU-
C LATED. NOTE THAT THERE ARE TWO POINTS AT THE INTERFACE.
C
      DIN1= 0.0
      NX=NXA+NXB
      THA=RC-THB
      R1=THB-RD
      THBI=THB
      THAI=THA
      DXB=THB/(NXB-1)
      IF(R1.GT.0.0)DXB=(THB-R1)/(NXB-2)
      DXA=THA/(NXA-1)
      TCONV=DMAX /RC**2
      XN(1) = 0.0
      XN(2)=DXB
      IF(R1.GT.0.0) XN(2)=R1
      DO 11  M=3,NXB
   11 XN(M) = XN(M-1) + DXB
      XN(NXB+1) =XN(NXB)
      NXAS=NXB+2
      DO 72 M=NXAS,NX
   72 XN(M)=XN(M-1)+DXA
      WRITE(6,114)     NXA,NXB,THB,THA,RC,RCF,RD,PCS
      NXB1=NXB+1
      NXB2=NXB+2
      NXB3=NXB+3
      NXBM1=NXB-1
      NXBM2=NXB-2
      NXBM3=NXB-3
      NXM1=NX-1
      NXO=NX
      NXAR=0
      NXBR=0
      XART=1.
      XBRT=1.
      NGRID=0
      WRITE(6,104)
      IF(NSP.GT.0) GO TO 411
      DO 481 I=1,NXBM1
  481 CSN(I  )=CPRIME
      CSN(NXB  )=CBA
      CSN(NXB1  )=CAB
      DO 41 I=NXB2,NX
      CSN(I  )=CZERO
   41 CONTINUE
      GO TO 410
C
C CONTROL IS TRANSFFERED TO THIS POINT FOR GRID CHANGE.
C
   12 DO 13 I=1,NX
   13 XSN(I)=XN(I)*RC
      NGRID=1
      WRITE(6,104)
      WRITE(6,107) TFINAL,THBF, THAF
      R1=R1*RC
      THB=THB*RC
      RD=THB-R1
      IF (CP(NX-2) .GT.CZERO1 ) RC=RC+RC-THB
      IF( RC.GT.RCF) RC=RCF
      THA=RC-THB
      IF(CP(4) .LT. CPRIME1) RD=RD*2
      R1=THB-RD
      IF(RD.GE.THB) R1=0.0
      IF(((CAB-CP(NX))/(CAB-CZERO)).LT.(0.3/ XART) .AND.(NXA.GT.5))
     1 NXAR=1
      IF(((CP(1)-CBA)/ (CPRIME-CBA)) .LT.(0.3/ XBRT).AND.(NXB.GT.5))
     1 NXBR=1
      IF(NXBR.EQ.0 .AND. NXAR.EQ.0)  GO TO 800
      NXO=NX
      IF(NXBR.EQ.1) NXB=NXB/2
      IF(NXAR.EQ.1) NXA=NXA/2
      IF(NXA.LT.5) NXA=5
      IF(NXB.LT.5) NXB=5
      NX=NXA+NXB
      NXB1=NXB+1
      NXB2=NXB+2
      NXB3=NXB+3
      NXBM1=NXB-1
      NXBM2=NXB-2
      NXBM3=NXB-3
      NXM1=NX-1
      NXAR=0
      NXBR=0
      XBRT= XBRT+2
      XART= XART+2
  800 CONTINUE
      DXB=(THB-R1)/(NXB-2)
      IF(R1.LT.DXB) R1=0.00
      IF( R1.LE.1.E-10 ) DXB= (THB-R1)/(NXB-1)
      DXA=THA/(NXA-1)
      TCONV=DMAX /RC**2
      XN(1)=0.0
      XN(2)=DXB
      IF(R1.GT.0.0) XN(2)=R1
      DO 14 M=3,NXB
   14 XN(M)=XN(M-1)+DXB
      XN(NXB1)=XN(NXB)
      DO 15 M=NXB2,NX
   15 XN(M)=XN(M-1)+DXA
  411 CSN(1  )=CP(1)
      DO 16 M=2,NXBM1
      I=1
   17 I=I+1
      IF(XN(M).GT.XSN(I) ) GO TO 17
      I1=I-1
   16 CSN(M)=CP(I1)-(CP(I1)-CP(I))/(XSN(I)-XSN(I1))*(XN(M)-XSN(I1))
      CSN(NXB)= CBA
      CSN(NXB1)=CAB
      DO 18 M=NXB2,NX
      IF(XN(M).GT.XSN(NXO)) GO TO 20
      I=1
   19 I=I+1
      IF(XN(M).GT.XSN(I)) GO TO 19
      I1=I-1
      CSN(M)=CP(I1)-(CP(I1)-CP(I))/(XSN(I)-XSN(I1)) *(XN(M)-XSN(I1))
      GO TO 18
   20 CSN(M  )=CP(NXO)
   18 CONTINUE
      IF( NGRID.GT.0 )  GO TO 624
  410 CONTINUE
      SUMI=0.
      CONST=1.
      IF(PCS.EQ. 1 ) CONST=  3.14159
      IF(PCS.EQ. 2 ) CONST= 4./3. * 3.14159
      DO 5    M=2,NX
      IF(M.EQ.NXB1)SUMIB=SUMI
    5 SUMI= SUMI+   0.5 * (CSN(M)+CSN(M-1))* CONST*(XN(M)**(PCS+1)-
     1XN(M-1)**(PCS+1))
  624 DO  625 M=1,NX
  625 CP(M)=CSN(M  )
      IF(NGRID.EQ.0) WRITE(6,132)
      IF(NGRID.EQ.1) WRITE(6,133)
      WRITE(6,115)
      WRITE(6,111) (I, XN(I),CSN(I), I=1,NX)
      DO  626 M=1,NX
  626 XN(M)=XN(M)/RC
      EOL=THB/RC
      THB=THB/RC
      THA=THA/RC
      DXA=DXA/RC
      DXB=DXB/RC
      R1=R1/RC
      IF(NGRID.GT.0) GO TO 23
      NP=NX
C
C NSOL      NUMBER OF DIFFUSION TIMES AT WHICH THE OUTPUT IS REQUIRED.
C ATIME     DIFFUSION TIME IN SECONDS AT WHICH OUTPUT IS REQUIRED.
C           THE PROGRAM RECYCLES BACK TO STATEMENT 38 TO READ A NEW DIFFUSION
C           TIME. NSOL DIFFUSION TIMES ARE READ IN ALTOGETHER.
C
      CPRIME1= CPRIME-0.0001
      CZERO1= CZERO +0.0001
      READ(5,101) NSOL
      LSOL = 0
   38 LSOL = LSOL+1
      KZ=0
      READ (5,112)  ATIME
   23 ITEST=0
   40 KZ=KZ+1
      CALL COEF
      DABAV= SQRT(DCAB*DMAX* DP(NXB2 ))
      DBAAV= SQRT(DCBA*DMAX* DP(NXBM1))
      DIDT= 1./ CDELTA*(DABAV/DMAX * (-CSN(NXB3  ) +4.* CSN(NXB2  )-3.
     1 *CAB)/2./DXA-DBAAV/DMAX * (CSN(NXBM2  )- 4.*CSN(NXBM1  ) +3.*CBA
     2   )/2./DXB)
      TIA=0.25*(DXA**2)  /DMAXA
      TIB=0.25*(DXB**2)  /DMAXB
      IF(CP(1) .LE. (CBA+1.E-06))TIB=1.E+70
      IF(CP(NX).GE. (CAB-1.E-06))TIA=1.E+70
      IF(TIA.GE.1.E+60 .AND. TIB.GE.1.E+60 )  GO TO 430
      TI= DMAX* AMIN1( TIA,TIB)
C
C DIN1,TOTAL DISTANCE THE INTERFACE HAS MOVED
C DINT,AMOUNT INTERFACE MOVED DURING TIME INCREMENT TI
C
      DINT= TI*DIDT
      IF(CP(1) .LE. (CBA+1.E-06))GO TO 47
      CP(1)=CSN(1  )+TI*GEF/2.0      /DXB**2*DP(1)*
     1(CSN(2  )-CSN(1  ))
C
C  DO LOOP 46 PERFORMS THE HEART OF THE CAL. FOR THE BETA PHASE
C
      DO 46 M=2,NXBM1
      C1= XN(M)/ EOL*DIDT*(CSN(M+1  )-CSN(M-1  )) /2. /DXB
      C2= DP(M) *(PCS/XN(M)*( CSN(M+1  )- CSN(M-1  ))/2./DXB +
     1 (CSN(M+1  ) -2.* CSN(M  ) +CSN(M-1  )) /DXB**2 )
      DJ= DP(M)*(ALOG(DP(M+1))- ALOG(DP(M-1)))
      C3= 0.25*        DJ         *  (CSN(M+1  )- CSN(M-1  ))  /DXB**2
      CP(M)= CSN(M    ) +TI* (C1+C2+C3)
      CSN(M-1  )=CP(M-1)
   46 IF(CP(M).GT.1.0) CP(M)=1.0
      CSN(NXBM1  )=CP(NXBM1)
C
C  DO LOOP 42 PERFORMS THE HEART OF THE CAL. FOR THE ALPHA-PHASE
C
   47 IF(CP(NX) .GE. (CAB-1.E-06))GO TO 48
      DO 42 M=NXB2,NXM1
      IF(CP(M-1).EQ.0.0) GO TO 421
      C1= (1. -XN(M))/ (1.-EOL)* (CSN(M+1  )-CSN(M-1  )) /2./DXA*DIDT
      C2= DP(M) *(PCS/XN(M)*( CSN(M+1  )- CSN(M-1  ))/2./DXA +
     1 (CSN(M+1  ) -2.* CSN(M  ) +CSN(M-1  )) /DXA**2 )
      DJ= DP(M)*(ALOG(DP(M+1))- ALOG(DP(M-1)))
      C3= 0.25*        DJ         *  (CSN(M+1  )- CSN(M-1  ))  /DXA**2
      CP(M)= CSN(M  ) +TI* (C1+C2+C3)
  421 IF(CP(M-1).LT.1.0E-05) CP(M)=0.0
   42 CSN(M-1  )=CP(M-1)
      CP(NX)= CSN(NX  ) +TI*2.*  DP(NX)*(CSN(NX-1  ) -CSN(NX  ))/DXA**2
      CSN(NXM1  )=CP(NXM1)
      CSN(NX  )=CP(NX)
   48 CONTINUE
      R1=R1+R1/THB*DINT
      THB=THB+DINT
      IF((THB*RC/THBI) .LT. 0.01)THB= 0.01*THBI /RC
      TI=TI/TCONV
      TFINAL=TFINAL+TI
      ITRIP=1
  802 CONTINUE
      THA=RC/RC-THB
      DXB=THB/NXBM1
      IF(R1.GT.0.0) DXB=(THB-R1)/NXBM2
      DXA=THA/(NXA-1)
      EOL=THB
      XN(1)=0.0
      XN(2)= DXB
      IF(R1.GT.0.0) XN(2)=R1
      DO 43 M=3,NXB
   43 XN(M)= XN(M-1) +DXB
      XN(NXB1)=XN(NXB)
      DO 44 M=NXB2,NX
   44 XN(M)= XN(M-1) + DXA
      IF(ITRIP.EQ.4)  GO TO 801
      SUMF=0.
      DO 60 M=2,NX
   60 SUMF= SUMF+0.5*(CSN(M)+CSN(M-1))* CONST*((XN(M)*RC) **(PCS+1)-
     1 (XN(M-1)*RC)**(PCS+1))
      DZI=  ((SUMI-SUMF)/CONST/CBA +(THB*RC)**(PCS+1 ))
      IF(DZI.LT.0.) DZI=0.
      THB= DZI**(1./ (1.+FLOAT(PCS)))/RC
      IF((THB*RC/THBI).LT.0.1)THB=.5**(1+PCS)*THB+ (1.-.5**(1+PCS))*EOL
      R1= R1+R1/EOL* (THB-EOL)
      ITRIP= ITRIP+1
      IF((THB*RC/THBI) .GE.0.10) ITRIP=4
      GO TO 802
  801 CONTINUE
      TNUM=1.0
      THAF=THA*RC
      THBF=THB*RC
      IF(R1.GT.0.0.AND.CP(4).LT.CPRIME1) GO TO 12
      IF(RC.LT.RCF.AND.CP(NX-2).GT.CZERO1 ) GO TO 12
      IF(((CP(1)-CBA)/ (CPRIME-CBA)) .LT.(0.3/ XBRT).AND.(NXB.GT.5))
     1 GO TO 12
      IF(((CAB-CP(NX))/(CAB-CZERO)).LT.(0.3/ XART) .AND.(NXA.GT.5))
     1 GO TO 12
      ITEST=ITEST+1
      IF (TFINAL.GE.ATIME) GO TO 430
      IF((THBF/THBI) .LT. 0.01)  GO TO 430
      GO TO 40
  430 CONTINUE
      TDIM= THBF-THBI
      DO 241 M=1,NX
  241 XN(M)=XN(M)*RC
      SUMF=0.
      DO 6    M=2,NX
      IF(M.EQ.NXB1)SUMFB=SUMF
    6 SUMF= SUMF+   0.5 * (CSN(M)+CSN(M-1))* CONST*(XN(M)**(PCS+1)-
     1XN(M-1)**(PCS+1))
      SUMBR= SUMFB/SUMIB
C
C THE FOLLOWING SEGMENT CONVERTS THE CONCENTRATIONS BACK TO ATOM FRACTION.
C
      DO 242 M=1,NX
  242 CP(M)=CP(M)*AMV/(BMV-CP(M)*(BMV-AMV))
      TFMIN=TFINAL/60.0
      CALL PLOTPF(CP,XN,NX,NA,XALM,XAL,YAL,DXN,DYN,XMIN,YMIN,XIN,YIN,
     1NCX,XCAP,NCY,YCAP,TYTLE,RESULT,LSOL,SPS)
      WRITE (6,104)
      WRITE (6,105)  ( ABLE(I), I=1,16 )
      WRITE(6,116) TSTART,TSMIN,     TFINAL,TFMIN,     THBI,THAI,
     1 THBF, THAF, RC, TDIM, SUMI, SUMF ,SUMBR
      WRITE(6,115)
      WRITE(6,111) (I, XN(I),CP (I), I=1,NX)
      IF(NO.NE.2) GO TO 58
      PUNCH 102,   (ABLE(I),I=1,16) , (FI(I),I=1,8 )
      PUNCH 103,   TSTART,TSMIN,     TFINAL,TFMIN,     THBI,THAI,
     1THBF,THAF,RC,DIN1         ,SUMI, SUMF
      PUNCH 101, NX
      PUNCH FI,   (CP(I),XN(I),I=1,NX)
   58 DO 260 I=1,NX
      XN(I)= XN(I) /RC
  260 CP(I)= CSN(I)
      IF((THBF/THBI) .LT. 0.01)  GO TO 992
      IF ( LSOL.NE.NSOL ) GO TO 38
  992 CONTINUE
      CALL NFRAME
      STOP
      END
      SUBROUTINE PLOTPF(CP,XN,NP,NA,XALM,XAL,YAL,DXN,DYN,XMIN,YMIN,
     1XIN,YIN,NCX,XCAP,NCY,YCAP,TYTLE,RESULT,LSOL,SPS)
C
C THIS SUBROUTINE IS WRITTEN TO PLOT THE CALCULATED CONCENTRATION PROFILES
C ON A CALCOMP PLOTTER WITH THE CDC COMPUTER SYSTEM AT NASA LANGLEY RESEARCH
C CENTER, HAMPTON, VA. TO RUN ON OTHER COMPUTING SYSTEMS, CHANGES IN THE PLOT
C STATEMENTS MAY BE REQUIRED.
C
      DIMENSION CP(1),XN(1),XCAP(1),YCAP(1),TYTLE(1),RESULT(1)
      DO 1 I=1,NP
      CP(I)=CP(I)*100.0
    1 XN(I)=XN(I)*1.0E4
      IF( NA.GT.1 .AND. LSOL.GT.1 )  GO TO 100
      TXIN=XIN
      WH=XAL/80.
      IF (WH.LT.0.08) WH=0.08
      IF (WH.GT.0.25) WH=0.25
      DXL=XAL/DXN
      DYL=YAL/DYN
      PA=3.0
    3 NX=XN(NP)/TXIN+0.95
      DXM=NX
      XDN=DXN
      IF(SPS.EQ.0.) XDN= AMAX1(DXN,DXM)
      XBL=DXL*XDN
      IF(XBL.GT.XALM) TXIN=PA/2.0*XIN
      PA=PA+1.0
      IF(SPS. EQ.0. AND. XBL.GT.XALM) GO TO 3
      NX=DXN+1.1
      NY=DYN+1.1
      TS=0.0
      TF=1.25*WH
      WK=1.4*WH
      YAL=AMIN1(YAL,9.0)
C
C          DRAW X-AXIS AND Y-AXIS
C
      CALL CALPLT(2.0,2.5,-3)
      CALL CALPLT(XBL,0.0,2)
      CALL CALPLT(XBL,YAL,2)
      CALL CALPLT(0.0,YAL,2)
      CALL CALPLT(0.0,0.0,2)
      VYP=-2.*WK
      CX=NCX
      CSP=XBL/2.-(CX/2.+1.)*WK
C
C          PLACE NUMBERS ALONG X-AXIS
C
      DO 5 I=1,NX
      RI=I-1
      TP=RI*DXL
      TV=XMIN+RI*TXIN+0.1
      IF (TV.GE.100.) VXP=TP-1.5*WK
      IF (TV.LT.100..AND.TV.GE.0.2) VXP=TP-WK
      IF (TV.LT.0.2) VXP=TP-0.5*WK
      CALL CALPLT(TP,TS,3)
      CALL CALPLT(TP,TF,2)
      IF(TV.LE.0.001) GO TO 5
      CALL NUMBER (VXP,VYP,WK,TV,0.,-1)
    5 CONTINUE
C
C          LABEL X-AXIS
C
      CALL SCRIBE(CSP,-7.*WH,1.4*WH,0.05,XCAP,0.,NCX,9)
      CY=NCY
      CSP=YAL/2.0-(CY/2.0+1.0)*WK
C
C          PLACE NUMBERS ALONG Y-AXIS
C
      DO 10 I=1,NY
      RI=I-1
      TP=RI*DYL
      VYP=TP-0.5*WK
      TV=YMIN+RI*YIN+0.1
      IF (TV.GE.100.) VXP=-4.*WK
      IF (TV.LT.100..AND.TV.GE.0.2) VXP=-3.*WK
      IF (TV.LT.0.2) VXP=-2.*WK
      CALL CALPLT(TS,TP,3)
      CALL CALPLT(TF,TP,2)
      IF(TV.LE.0.0) GO TO 10
      CALL NUMBER (VXP,VYP,WK,TV,0.,-1)
   10 CONTINUE
      XAL2=XBL/4.0
      XAL3=XAL2+2.0
C
C          LABEL Y-AXIS
C
      CALL SCRIBE(-7.*WH,CSP,1.4*WH,0.05,YCAP,90.,NCY,9)
      CALL SCRIBE(XAL2,4.75,1.40*WH,0.05,RESULT,0.,20,9)
      CALL SCRIBE(XAL2,4.50,1.40*WH,0.05,TYTLE,0.0,80,9)
  100 CONTINUE
      DO 15 I=1,NP
      X =( XN(I)-XMIN)*XBL/(TXIN*DXN)
      Y=( CP(I  )-YMIN)*YAL/(YIN*DYN)
      IF (Y.LT.0.) Y=0.
      IF (Y.GT.YAL) Y=YAL
      IF(X.GT.XBL) X=XBL
      IF (X.LE.0..OR.Y.LE.0.) CALL CALPLT(X,Y,3)
      IF(X.GT.0.0.AND.Y.GT.0.0) CALL CALPLT(X,Y,2)
   15 CONTINUE
      DO 2 I=1,NP
      CP(I)=CP(I)/100.0
    2 XN(I)=XN(I)/1.0E4
      CALL CALPLT(0.0,0.0,3)
      IF(NA.EQ.1)  CALL NFRAME
      RETURN
      END
      SUBROUTINE COEF
      COMMON/STORE/NX,CN(51 ),DN(51), ND,CP(501),DP(501) ,DMAX
C
C DIFFUSION COEFFICIENT CALCULATION BY LOGERTHIMIC INTERPOLATION.
C
      DO 4 M=1,NX
      DO 2 I=2,ND
      IF(CP(M).LT.CN(I)) GO TO 2
      DP(M)=DN(I)*(DN(I-1)/DN(I))**((CP(M)-CN(I))/(CN(I-1)-CN(I)))
      GO TO 4
    2 CONTINUE
    4 DP(M)= DP(M)/DMAX
      RETURN
      END
